Coronavirus and paramyxovirus shedding by bats in a cave and buildings in Ethiopia
- Setting up the workspace and loading data
- Sampling events
- Viral families tested
- Species of bats tested
- Number of bats testing positive
- Figure 1. Sampling sites
- Table 1. Summary of bats sampled and viral testing results
- Figure 2. Summary by sampling month, season, bat species, and sampling interface/location
- Figure 3. Viral findings by bat species and sampling interface
- Viruses detected
- Coinfections
- Virus detection by specimen type
- Table 2. Viral prevalence in bats by sex, age class, season, species, and specimen type
- Table 3. Association between bats testing positive for coronaviruses or paramyxoviruses and their species, sex, age, and season in which they were sampled
Code for manuscript.
Setting up the workspace and loading data
# Packages
# ==============================================================================
library(tidyr)
library(ggplot2)
library(knitr)
library(DT)
library(stringr)
library(xtable)
library(hrbrthemes)
library(dplyr)
library(forcats)
library(finalfit)
library(dplyr)
library(stringr)
library(janitor)
library(extrafont)
library(cowplot)
library(knitr)
library(parameters)
library(here)
# Load fonts
extrafont::loadfonts(quiet = TRUE)
# Load data
# ==============================================================================
bats_res = readRDS("data/Ethiopia_bat-info_test-summary.RDS")
bats_res$SiteName = trimws(gsub("-\\<Bat\\>|\\<Bat\\>", "", bats_res$SiteName))
bats_res$ConcurrentSamplingSite = bats_res$SiteName
tests = readRDS("data/Ethiopia_bat-test-info.RDS")
tests$SiteName = trimws(gsub("-\\<Bat\\>|\\<Bat\\>", "", tests$SiteName))
tests$Interface = ifelse(tests$ScientificName %in% "Rhinopoma hardwickii",
"Cave", "Building")
events = readRDS("data/Ethiopia_bat-sampling-events.RDS")
df = bats_res
df$Family = str_to_title(df$Family)
df$AgeClass = word(df$AgeClass, 1)
df$Interface = ifelse(df$ScientificName %in% "Rhinopoma hardwickii",
"Cave", "Building")
date_order = unique(format(as.Date(sort(df$EventDate), format = "%Y-%m-%d"),
"%Y-%B"))
df$EventMY = format(as.Date(df$EventDate, format = "%Y-%m-%d"), "%Y-%B")
df$EventMY = factor(df$EventMY, levels = date_order, ordered = TRUE)
tests = left_join(tests, df[ , c("AnimalID", "Family")], by = "AnimalID")Sampling events
df %>%
group_by(SiteName, EventDate, SeasonModelled) %>%
tally(name = "No. of bats sampled") %>%
group_by(EventDate) %>%
arrange(SiteName, EventDate) %>%
adorn_totals("row") %>%
datatable(options = list(dom = 't'), rownames = FALSE)Species of bats tested
df %>%
select(Family, ScientificName) %>%
group_by(Family, ScientificName) %>%
tally() %>%
adorn_totals("row") %>%
datatable(options = list(dom = 't'), rownames = FALSE)Number of bats testing positive
Number of bats tested = 589
Number of bats testing positive = 99
Percent bats testing positive = 16.8081494
Figure 1. Sampling sites
Made with ArcGIS 10.1
out_file = file.path(here(), "manuscript_figures")
myimages = list.files(out_file, pattern = glob2rx("Figure-1*png"),
full.names = TRUE)
include_graphics(myimages)
Figure 2. Summary by sampling month, season, bat species, and sampling interface/location
Figure subcomponent 1: positives by sampling month and season
p_df = df %>%
group_by(EventMY, is_positive) %>%
arrange(EventMY) %>%
summarize(n = n()) %>%
mutate(is_positive = ifelse(is_positive, "Virus Positive", "Virus Negative")) %>%
mutate(is_positive = factor(is_positive, levels = c("Virus Positive", "Virus Negative")))
levels(p_df$EventMY) = gsub("-", "\n", levels(p_df$EventMY))
event_plot = ggplot(p_df, aes(x = EventMY, y = n, fill = is_positive, label = n)) +
geom_bar(stat="identity", width = 0.5) +
scale_fill_manual(values = c("#BD4F2C", "#2C42BD"), name = "Test result") +
theme_ipsum(base_size = 12) +
ylab("Number of bats") +
xlab("") +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
legend.position = "top") +
scale_y_continuous(expand = c(0,1))
p_season_df = df %>%
group_by(EventMY, SeasonModelled) %>%
arrange(EventMY) %>%
summarize(n = n())
levels(p_season_df$EventMY) = gsub("-", "\n", levels(p_season_df$EventMY))
event_plot = event_plot + geom_text(data = p_season_df,
aes(label = SeasonModelled,
x = EventMY, y = n,
fill = NULL),
nudge_y = 2,
size = 3)
cowplot::save_plot("results/figs/event-plot.pdf", event_plot, base_height = 6)
cowplot::save_plot("results/figs/event-plot.png", event_plot, base_height = 6)
event_plot
Figure subcomponent 2: bat species tested
bat_plot = df %>%
select(ScientificName, is_positive) %>%
group_by(ScientificName, is_positive) %>%
summarize(n = n()) %>%
arrange(desc(n)) %>%
as.data.frame %>%
mutate(ScientificName = factor(ScientificName,
levels = c(
"Neoromicia cf. somalica",
"Mops midas",
"Chaerephon pumilus",
"Rhinopoma hardwickii"
))) %>%
ggplot( aes(x = ScientificName, y= -n,
fill = as.factor(-is_positive), width = 0.75)) +
geom_bar(stat="identity") +
scale_fill_manual(values = c("#BD4F2C", "#2C42BD"), name = "") +
scale_x_discrete(name = "", position = "top") +
scale_y_continuous(name = "No. of bats",
breaks = seq(0, -300, by = -50), # y axis values (before coord_flip)
labels = seq(0, 300, by = 50)) +
coord_flip() +
theme_ipsum(base_size = 12, plot_title_size = 12) +
theme(
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
legend.position="none") +
xlab("") +
ylab("Number of bats") +
labs(title = "Species tested")
cowplot::save_plot("results/figs/bat-plot.pdf", bat_plot, base_height = 3)
bat_plot
Figure subcomponent 3: Pie charts for species-event timeline
sp_df = df %>%
select(ScientificName, is_positive, EventDate) %>%
group_by(ScientificName, is_positive, EventDate) %>%
summarize(n = n()) %>%
arrange(desc(n)) %>%
mutate(sp_e = paste0(EventDate, "_", ScientificName)) %>%
as.data.frame
sp_l = split(sp_df, sp_df$sp_e)
# ref:
# http://www.sthda.com/english/wiki/ggplot2-pie-chart-quick-start-guide-r-software-and-data-visualization
blank_theme = theme_minimal()+
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
panel.border = element_blank(),
panel.grid=element_blank(),
axis.ticks = element_blank(),
plot.title=element_text(size=14, face="bold"),
legend.position = "none"
)
out_dir = "results/figs/event_pie"
if( !dir.exists(out_dir)) dir.create(out_dir)
invisible(lapply(sp_l, function(sp, out_dir) {
if(nrow(sp) > 1) {
sp_p = sp %>%
ggplot(aes(x = "", y = n, fill = as.factor(-is_positive))) +
geom_bar(width = 1, stat="identity") +
scale_fill_manual(values = c("#BD4F2C", "#2C42BD"), name = "") +
coord_polar("y", start=0) +
theme_minimal() +
blank_theme
} else {
sp_p = sp %>%
ggplot(aes(x = "", y = n, fill = "#B9B7B6")) +
geom_bar(width = 1, stat="identity") +
scale_fill_manual(values = "#2C42BD", name = "") +
coord_polar("y", start=0) +
theme_minimal() +
blank_theme
}
out_name_pdf = file.path(out_dir, paste0(gsub(" ", "_", unique(sp$sp_e)), ".pdf"))
out_name_png = file.path(out_dir, paste0(gsub(" ", "_", unique(sp$sp_e)), ".png"))
cowplot::save_plot(out_name_pdf, sp_p, base_height = .5)
cowplot::save_plot(out_name_png, sp_p, base_height = .5)
}, out_dir))out_dir = file.path(here(), "results/figs/event_pie")
myimages = list.files(out_dir, pattern = ".png", full.names = TRUE)
include_graphics(myimages)
Final Figure 2 compiled with above subcomponents in Affinity Designer 1.8.3
out_file = file.path(here(), "manuscript_figures")
myimages = list.files(out_file, pattern = glob2rx("Figure-2*png"),
full.names = TRUE)
include_graphics(myimages)
Figure 3. Viral findings by bat species and sampling interface
Figure subcomponent 1: bar chart
df_p =
df %>%
mutate(SeasonModelled = factor(SeasonModelled),
SeasonModelled = fct_rev(SeasonModelled),
ScientificName = ifelse(ScientificName %in% "Rhinopoma hardwickii",
paste0(ScientificName, " (Cave)"),
paste0(ScientificName, " (Building)")))
virus_plot =
df_p %>%
select(pos_virus, pos_virus_n, pos_virus_family,
ScientificName, SeasonModelled) %>%
drop_na() %>%
mutate(pos_virus = gsub("; ", "\nand ", pos_virus)) %>%
group_by(pos_virus, pos_virus_n, pos_virus_family,
ScientificName, SeasonModelled) %>%
summarize(n = n()) %>%
arrange(pos_virus_family, desc(n)) %>%
as.data.frame %>%
mutate(pos_virus = factor(pos_virus, unique(pos_virus)),
pos_virus = fct_rev(pos_virus)) %>%
ggplot(aes(x = reorder(pos_virus, n), y=n, fill = ScientificName, label = n) ) +
facet_wrap( ~ SeasonModelled ) +
scale_fill_manual(values = c("#CC6677", "#44AA99", "#DDCC77", "#AA4499"),
name = "Bat species and\nsampling interface") +
geom_bar(stat="identity") +
ylim(0, 80) +
# geom_text( size = 3, nudge_y = 4) +
coord_flip() +
theme_ipsum(base_size = 10) +
theme(
axis.text.y = element_text(size = 8),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
legend.position="right",
legend.text = element_text(face = "italic")
) +
xlab("") +
ylab("Number of bats sampled") +
labs(title = "Viruses detected in bats")
cowplot::save_plot("results/figs/virus-plot.pdf",
virus_plot, base_height = 6)
virus_plot
Figure subcomponent 2: Extra legend information
df %>%
group_by(ScientificName, SeasonModelled, is_positive) %>%
tally() %>%
mutate(is_positive = ifelse(is_positive, "Positive", "Negative")) %>%
pivot_wider(names_from = is_positive, values_from = n) %>%
adorn_totals("col") %>%
mutate_all(~replace_na(., 0)) %>%
mutate(pos_percent = round((Positive/Total * 100), 2)) %>%
mutate(info = paste0(pos_percent, "% (", Positive, "/", Total, ")")) %>%
select(ScientificName, SeasonModelled, info) %>%
pivot_wider(names_from = SeasonModelled, values_from = info) %>%
datatable(options = list(dom = 't'))Final Figure 3 compiled with above subcomponents in Affinity Designer 1.8.3
out_file = file.path(here(), "manuscript_figures")
myimages = list.files(out_file, pattern = glob2rx("Figure-3*png"),
full.names = TRUE)
include_graphics(myimages)
Viruses detected
tbl = xtable(table(unlist(strsplit(df$pos_virus[ df$is_positive ], "; "))))
names(tbl) = "No. of bats testing positive for virus"
datatable(tbl)Coinfections
df %>%
filter( is_positive & pos_virus_n > 1 ) %>%
select(SiteName, ScientificName, pos_virus) %>%
group_by(SiteName, ScientificName, pos_virus) %>%
summarize(n_bats = n()) %>%
datatable(colnames = c("SiteName", "ScientificName", "Virus coinfection",
"No. of bats"))Virus detection by specimen type
Breakdown of positive specimen types (at individual bat level)
Rectal swabs were the most common specimen type to yield a positive virus sequence (92/99 bats), including 27 bats from whom sequences were confirmed in both rectal and oral swabs. In addition, seven bats had viruses detected only via their oral swabs.
df %>%
filter(is_positive) %>%
group_by(pos_specimen_type) %>%
tally() %>%
adorn_totals("row") %>%
datatable(options = list(dom = 't'))Positive specimen types for bats with > 1 virus detected
ids = df$AnimalID[ which(df$pos_virus_n > 1) ]
multi_bats = unique(
tests[ tests$AnimalID %in% ids & !is.na(tests$VirusGroup),
c("AnimalID", "SpecimenType", "VirusGroup")]
)
tbls = vector("list", 3)
for(i in seq_along(ids)) {
tbls[[ i ]] =
multi_bats %>%
filter(AnimalID %in% ids[ i ]) %>%
datatable(options = list(dom = 't'))
}
tbls[[1]]tbls[[2]]tbls[[3]]Virus group by specimen type (at individual bat level)
df %>%
filter(is_positive) %>%
select(pos_virus, pos_specimen_type) %>%
group_by(pos_virus, pos_specimen_type) %>%
tally() %>%
as_data_frame() %>%
spread(pos_specimen_type, n) %>%
mutate_all(~replace_na(., 0)) %>%
adorn_totals(c("row", "col")) %>%
datatable(options = list(dom = 't'),
colnames = c("Virus", "Oral swab only", "Oral and rectal swab",
"Rectal swab only", "Total"),
rownames = FALSE)